
//This is the meat-and-potatoes of the program. Here the actual
//simulation code is defined, drawing from two abstract classes
//mode_method and thread_method, which allow the simulator object
//to run the same code with different settings.

#ifndef SIMULATOR_H
#define SIMULATOR_H

#include "grid.hpp"
#include "settings.hpp"
#include "bitmap.hpp"
#include <vector>
#include <memory>
#ifdef CAVITY_MULTITHREADING_AVAILABLE
  #include <boost/thread/barrier.hpp>
  #include <boost/thread/mutex.hpp>
  #include <boost/thread/thread.hpp>
  #include "sync_func.hpp"
#endif

///////////////Interface Types, to avoid extensive switch statements///////////////
//////////////////////////////////////////////////////////////////////////////////

class thread_method //multi-thread or single-thread
{
  public:
    virtual bool wait() = 0;
    virtual grid::data_type get_sum(const grid::data_type & piece) = 0;
    virtual grid::data_type get_max(const grid::data_type & piece) = 0;
    virtual grid::data_type get_min(const grid::data_type & piece) = 0;
    virtual grid::region & region(std::vector<grid::region>::size_type index = 0) = 0;
    virtual grid::boundary & boundary(std::vector<grid::boundary>::size_type index = 0) = 0;
};

class mode_method //one mode or many modes
{
  public:
    mode_method(thread_method & tmeth_) : tmeth(tmeth_) {}
    virtual void subtract(grid::region & reg) = 0;
    virtual void adjust_boundary(const grid::boundary & bound, grid & gref) = 0;
    virtual void update_sumofsq(grid & g,grid::region & reg) = 0;
			//Calculate and store the sum of squares of the just-finished mode;
            //so that when (if) it is subtracted from higher modes later, its sum of squares need not
            //be calculated again. If we're dealing with only one mode ( e.g. typeid(*mmeth) == typeid(low_mode()) ),
            //then this function will do nothing.
  protected:
    thread_method & tmeth;
};

class mode_type //TE or TM
{
    public:
      mode_type(thread_method & tmeth_) : tmethod(tmeth_) {}
      virtual void adjust_boundary(const grid::boundary & bound, grid & gref) = 0;
    protected:
      thread_method & tmethod;
};

////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////

class TE_mode : public mode_type
{
    public:
      TE_mode(thread_method & tmeth_) : mode_type(tmeth_) {}
      void adjust_boundary(const grid::boundary & bound, grid & gref);
};

class TM_mode : public mode_type
{
    public:
        TM_mode(thread_method & tmeth_) : mode_type(tmeth_) {}
        void adjust_boundary(const grid::boundary & bound, grid & gref) { /*empty*/ }
};

class high_mode : public mode_method
{
    vec_of_grids & modes;
    std::vector<grid::data_type> sosq_vec; //store the sums of squares
  public:
    high_mode(thread_method & tmeth_, vec_of_grids & modes_)
        : mode_method(tmeth_), modes(modes_) { /*empty*/ }
    void subtract(grid::region & reg); //uses the modes and sums of squares
    void update_sumofsq(grid & g, grid::region & reg);
        //associate with each mode it's sum of squares, after we're done calculating the mode.
};

struct low_mode : public mode_method
{
    low_mode(thread_method & tmeth_) : mode_method(tmeth_) { /*empty*/ }
    void subtract(grid::region & reg) { /*empty*/ }
    void update_sumofsq(grid & g,grid::region & reg) { /*empty*/ }
};

struct high_TM : public high_mode, public TM_mode
{
    high_TM(thread_method & tmeth_, vec_of_grids & modes_) : high_mode(tmeth_, modes_), TM_mode(tmeth_) {}
    void adjust_boundary(const grid::boundary & bound, grid & gref) { TM_mode::adjust_boundary(bound,gref); }
};

struct low_TM : public low_mode, public TM_mode
{
    low_TM(thread_method & tmeth_) : low_mode(tmeth_), TM_mode(tmeth_) {}
    void adjust_boundary(const grid::boundary & bound, grid & gref) { TM_mode::adjust_boundary(bound,gref); }
};

class high_TE : public high_mode, public TE_mode
{
      void partition_boundary(const grid::boundary & b, grid::size_type nthreads);
    public:
      high_TE(thread_method & tmeth_, vec_of_grids & modes_) : high_mode(tmeth_, modes_), TE_mode(tmeth_) {}
      //grid::boundary & get_bound(std::vector<grid::boundary>::size_type index) { return bounds[index]; }
      void adjust_boundary(const grid::boundary & bound, grid & gref) { TE_mode::adjust_boundary(bound,gref); }
};

struct low_TE : public low_mode, public TE_mode
{
    low_TE(thread_method & tmeth_) : low_mode(tmeth_), TE_mode(tmeth_) {}
    void adjust_boundary(const grid::boundary & bound, grid & gref) { TE_mode::adjust_boundary(bound,gref); }
};

class single_thread : public thread_method
{
    public:
        bool wait() { return true; }
        grid::data_type get_sum(const grid::data_type & piece) { return piece; }
        grid::data_type get_max(const grid::data_type & piece) { return piece; }
        grid::data_type get_min(const grid::data_type & piece) { return piece; }
        grid::region & region(std::vector<grid::region>::size_type index)
        {
            if(index) throw std::string("single_thread requested region index other than 0");
            return r;
        }
        grid::boundary & boundary(std::vector<grid::boundary>::size_type index)
        {
            throw std::string("In single_thread::boundary()... I don't think this function should be called!");
            return dummy;
        }
    private:
        grid::region r;
        grid::boundary dummy;
};

#ifdef CAVITY_MULTITHREADING_AVAILABLE
class multi_thread : public thread_method
{
        boost::barrier bar;
        sync_func< grid::data_type, thread_max<grid::data_type> > sync_max;
        sync_func< grid::data_type, thread_min<grid::data_type> > sync_min;
        sync_func< grid::data_type, thread_sum<grid::data_type> > sync_sum;
        std::vector<grid::region> r;
        std::vector<grid::boundary> b;
        template<class vec_T>
        void partition(const vec_T & in, std::vector<vec_T> & out, grid::size_type nthreads);
    public:
        multi_thread(unsigned nthreads, grid::region & reg, grid::boundary & bound)
            : bar(nthreads), sync_max( nthreads, std::numeric_limits<grid::data_type>::min() ),
              sync_min( nthreads, std::numeric_limits<grid::data_type>::max() ), sync_sum( nthreads, 0.0 )
            {
                partition(reg, r, nthreads);
                partition(bound, b, nthreads);
            }
        bool wait() { return bar.wait(); }
        grid::data_type get_sum(const grid::data_type & piece) { return sync_sum(piece); }
        grid::data_type get_max(const grid::data_type & piece) { return sync_max(piece); }
        grid::data_type get_min(const grid::data_type & piece) { return sync_min(piece); }
        grid::region & region(std::vector<grid::region>::size_type index) { return r.at(index); }
        grid::boundary & boundary(std::vector<grid::boundary>::size_type index) { return b.at(index); }
};
#endif

#ifdef CAVITY_MULTITHREADING_AVAILABLE
template<class vec_T> //vec_T is an object that has interface like std::vector.
void multi_thread::partition(const vec_T & in, std::vector<vec_T> & out, grid::size_type nthreads)
{
    std::vector<vec_T> & partition_v = out;
    partition_v.clear();
    uword insz = num_points(in);
	uword fmt_old = 0; //'fmt' for "first more than"
	uword fmt_new = 0;
	for(unsigned int i = 0; i < nthreads; ++i)
	{
		fmt_new = first_more_than( in, insz/nthreads, fmt_old+1 );
		partition_v.push_back( vec_T(in.begin()+fmt_old, in.begin()+fmt_new) );
		fmt_old = fmt_new;
	}
}
#endif
class simulator
{
    public:
        simulator(thread_method & tmeth_, mode_method & mmeth_, grid & g_,
                  grid & lap_, grid::region & reg_, grid::boundary & bound_,
                  uword n_, grid::data_type dt_, uword sub_period_ = settings::sub_period,
                  uword norm_period_ = settings::norm_period);
        void operator()();
        void low_TM_single_thread();
    private:
        void normalize(grid::data_type maximum = 1.0, grid::data_type scale = 1.0);
        thread_method & tmeth;	//wait(), get_sum(), get_max(), get_min()
        mode_method & mmeth;	//subtract(), update_sumofsq()
        grid::region & reg;
        grid::boundary & bound;
        grid & g;
        functions::calc_lap lapf;
        functions::relax relaxf;
        uword n;
        uword subtraction_period;
        uword normalization_period;
};

#endif //SIMULATOR_H

